/* 
Scenic Graphics Library
Copyright (C) 2007 Jouni Tulkki

This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.

This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USApackage scenic;
*/

#include <math.h>
#include <stdlib.h>


#define INTEGRAL_SIZE		32


struct LookupTableContext
{
	double (*func)(double, double, void*);
	void *data;
	double integral[INTEGRAL_SIZE][INTEGRAL_SIZE];
	double integral_2[INTEGRAL_SIZE][INTEGRAL_SIZE];
};

void calc_integral(double (*func)(double, double, void*), void *data, double *integral, 
				   double *integral_2, int w, int h)
{
	int x, y;
	double sum;

	for(y = 0; y < h; y++)
	{
		sum = 0;
		for(x = 0; x < w; x++)
		{
			sum += func((x + 0.5) / w, (y + 0.5) / (double)h, data);
			integral_2[x + y * w] = sum / w;
		}
	}

	for(x = 0; x < w; x++)
	{
		sum = 0;
		for(y = 0; y < h; y++)
		{
			sum += integral_2[x + y * w];
			integral[x + y * w] = sum / h;
		}
	}
}

static double step_function(double x)
{
	if(x >= 0.0)
		return 1.0;
	return 0.0;
}

static double calc_integral_2y(LookupTableContext *ctx, double y, int x)
{
	int iy = (int)floor(y);
	double ay = y - floor(y);

	if(y < 0)
		return ctx->integral_2[0][x];
	if(y >= INTEGRAL_SIZE - 1)
		return ctx->integral_2[INTEGRAL_SIZE - 1][x];

	return ctx->integral_2[iy][x] * (1.0 - ay) + ctx->integral_2[iy + 1][x] * ay;
}

static double calc_integral_2(LookupTableContext *ctx, double y, double x)
{
	int ix = (int)floor(x);
	double ax = x - floor(x);
	double offset = step_function(x - INTEGRAL_SIZE / 2) * 
		calc_integral_2y(ctx, y, INTEGRAL_SIZE - 1);

	if(x < 0 || x >= INTEGRAL_SIZE)
		return 0;
	if(ix == 0)
		return offset - calc_integral_2y(ctx, y, 0) * ax;
	return offset - (calc_integral_2y(ctx, y, ix - 1) * (1.0 - ax) + 
		calc_integral_2y(ctx, y, ix) * ax);
}

static double calc_integral(LookupTableContext *ctx, double y)
{
	int iy = (int)floor(y);
	double ay = y - floor(y);
	double offset = step_function(y - INTEGRAL_SIZE / 2);

	if(y < 0 || y > INTEGRAL_SIZE)
		return 0;
	if(iy == 0)
		return offset - ctx->integral[0][INTEGRAL_SIZE - 1] * ay;
	return offset - (ctx->integral[iy - 1][INTEGRAL_SIZE - 1] * (1.0 - ay) + 
		ctx->integral[iy][INTEGRAL_SIZE - 1] * ay);
}
#include <stdio.h>

static double calc_value_for_point(LookupTableContext *ctx, 
								   double angle, double x, double y)
{
	double sina = sin(angle);
	double cosa = cos(angle);
	double y2;
	double sum = 0.0;
	double ry, rx;

	if(sina > 0.0001)
	{
		for(y2 = 0.125; y2 < INTEGRAL_SIZE; y2 += 0.25)
		{
			ry = y + y2 / (double)INTEGRAL_SIZE - 0.5;
			rx = ry * (cosa / sina) - x;

			if(ry >= 0.0)
			{
				sum -= calc_integral_2(ctx, y2 - 0.5, (rx + 0.5) * INTEGRAL_SIZE);
			}
		}
		sum /= INTEGRAL_SIZE * 4;
	}
	if(cosa != 0.0 && x * cosa > 0.0)
	{
		ry = x * (sina / cosa) - y;
//		if(y2 >= 0 && y2 < INTEGRAL_SIZE)
		{
			sum += calc_integral(ctx, (ry + 0.5) * INTEGRAL_SIZE);
//			sum += step_function(y2 - INTEGRAL_SIZE / 2) - ctx->integral[y2][INTEGRAL_SIZE - 1];
		}
	}

	return sum;
}

static double calc_inverse_value_for_point(LookupTableContext *ctx,
										   double angle, double x, double y)
{
	double sina = sin(angle);
	double cosa = cos(angle);

	return calc_value_for_point(ctx, angle, x + cosa * 10, y + sina * 10) -
		calc_value_for_point(ctx, angle, x, y);
}

static void calc_lookup_table_for_angle(LookupTableContext *ctx, 
						 double *lookup_table, int size, double angle)
{
	int s = size - 3;
	double sina = sin(angle);
	double cosa = cos(angle);
	double y, x;
	int i, j;
	double m = 1.0 / (s - 1.0);

	printf("Angle: %lf\n", angle);
	for(i = 0; i < s + 1; i++)
	{
		y = i * m - 0.5;
		for(j = 0; j < s + 1 - i; j++)
		{
			x = j * m - 0.5;
			lookup_table[i * size + j] = (calc_value_for_point(ctx, angle, x - 0.001, y) + calc_value_for_point(ctx, angle, x - 0.001, y)) / 2.0;
		}
	}
	for(i = 0; i < s + 1; i++)
	{
		y = i * m - 0.5;
		for(j = 0; j < s + 1 - i; j++)
		{
			x = j * m - 0.5;
			lookup_table[(size - 1 - i) * size + (size - 1 - j)] = (calc_value_for_point(ctx, angle, x + 0.001, y) + calc_value_for_point(ctx, angle, x + 0.001, y)) / 2.0;
		}
	}
	for(i = 0; i < size; i++)
	{
		for(j = 0; j < size; j++)
		{
			if(abs(i + j - (size - 1)) <= 1)
			{
				if(i > j)
				{
					y = (i + 0.0001 - 3) * m - 0.5 + sina * 10;
					x = (j - 0.0001 - 0) * m - 0.5 + cosa * 10;
					lookup_table[i * size + j] = calc_value_for_point(ctx, angle, x, y);
					if(y * cosa - x * sina < 0.0)
						lookup_table[i * size + j] -= 1;
				}
				else if(i < j)
				{
					y = (i - 0.0001 - 0) * m - 0.5 + sina * 10;
					x = (j + 0.0001 - 3) * m - 0.5 + cosa * 10;
					lookup_table[i * size + j] = calc_value_for_point(ctx, angle, x, y);
					if(y * cosa - x * sina > 0.0)
						lookup_table[i * size + j] += 1;
				}
			}
		}
	}
}

void calcLookupTable(double (*func)(double, double, void*), void *data, 
					 double *lookup_table, int size, int d)
{
	LookupTableContext *ctx = (LookupTableContext*)malloc(sizeof(LookupTableContext));
	int i;
	double a, angle;

	calc_integral(func, data, &ctx->integral[0][0], &ctx->integral_2[0][0], 
		INTEGRAL_SIZE, INTEGRAL_SIZE);
	ctx->func = func;
	ctx->data = data;

	for(i = 0; i < d; i++)
	{
		a = i / (d - 1.0);
		if(a < 0.5)
			angle = atan2(a * 2, 1);
		else
			angle = atan2(1, (1.0 - a) * 2);
		calc_lookup_table_for_angle(ctx, lookup_table + i * size * size, size, angle);
	}

	free(ctx);
}

static void calc_raw_lookup_table_for_angle(LookupTableContext *ctx, 
						                    double *lookup_table, int size, double angle)
{
	int s = size;
	double y, x;
	int i, j;

	for(i = 0; i < size; i++)
	{
		for(j = 0; j < size; j++)
		{
			y = i / (s - 1.0) - 0.5;
			x = j / (s - 1.0) - 0.5;

			lookup_table[i * size + j] = (calc_value_for_point(ctx, angle, x + 0.0001, y + 0.0001) +
				calc_value_for_point(ctx, angle, x - 0.0001, y - 0.0001)) / 2;
		}
	}
}


void calcRawLookupTable(double (*func)(double, double, void*), void *data, 
					 double *lookup_table, int size, int d)
{
	LookupTableContext *ctx = (LookupTableContext*)malloc(sizeof(LookupTableContext));
	int i;
	double a, angle;

	calc_integral(func, data, &ctx->integral[0][0], &ctx->integral_2[0][0], 
		INTEGRAL_SIZE, INTEGRAL_SIZE);
	ctx->func = func;
	ctx->data = data;

	for(i = 0; i < d; i++)
	{
		a = i / (d - 1.0);
		if(a < 0.5)
			angle = atan2(a * 2, 1);
		else
			angle = atan2(1, (1.0 - a) * 2);
		calc_raw_lookup_table_for_angle(ctx, lookup_table + i * size * size, size, angle);
	}

	free(ctx);
}

